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This paper presents a CFD model for simulating the self-pressurization of a large scale 
liquid hydrogen storage tank. In this model, the kinetics-based Schrage equation is used to 
account for the evaporative and condensing interfacial mass flows. Laminar and turbulent 
approaches to modeling natural convection in the tank and heat and mass transfer at the 
interface are compared. The flow, temperature, and interfacial mass fluxes predicted by 
these two approaches during tank self-pressurization are compared against each other. The 
ullage pressure and vapor temperature evolutions are also compared against experimental 
data obtained from the MHTB self-pressurization experiment. A CFD model for cooling 
cryogenic storage tanks by spraying cold liquid in the ullage is also presented. The Euler- 
Lagrange approach is utilized for tracking the spray droplets and for modeling interaction 
between the droplets and the continuous phase (ullage). The spray model is coupled with the 
VOF model by performing particle tracking in the ullage, removing particles from the ullage 
when they reach the interface, and then adding their contributions to the liquid. Droplet- 
ullage heat and mass transfer are modeled. The flow, temperature, and interfacial mass flux 
predicted by the model are presented. The ullage pressure is compared with experimental 
data obtained from the MHTB spray bar mixing experiment. The results of the models with 
only droplet/ullage heat transfer and with heat and mass transfer between the droplets and 
ullage are compared. 
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I. Introduction 

Efficient cryogenic fluid management of high energy propellants are in the critical path of nearly all the 
NASA future human exploration mission scenarios 1 . A key technology challenge is long duration storage of 
cryogens in space, where heat leaks through the insulation and various conduction paths cause tank self- 
pressurization. For short term operations, thrusters have traditionally been used to settle the propellant and relieve 
the tank pressure by venting. However, the added propellant and hardware weight (cost) to accommodate the 
increasing number of settling/venting sequences that becomes necessary in long duration missions will be quite 
prohibitive. Thus, maintaining tank pressure while minimizing the boil-off loss of propellants through an active 
pressure control mechanism has become a significant challenge associated with long term storage of cryogens in 
microgravity". 

Among various pressure control strategies that have been proposed and tested, the spray-bar 
thermodynamic vent system (TVS) has emerged as a promising mechanism that enables tank pressure control 
through venting without resettling. The key components of a TVS include a Joule-Thompson (J-T) expansion, a 
two-phase heat exchanger, and a spray -bar mixing pump system which, working together, enable destratification and 
extraction of heat from the tank with minimized liquid cryogen losses 3 during pressure control. Unfortunately, 
implementation and optimization of the spray-bar TVS for future mission applications has been slow due to a lack of 
opportunities for on-orbit testing and technology demonstration. In order to partially alleviate this situation, it was 
decided to use comprehensive ground testing together with computational model development to aid the storage tank 
design. In this scheme, the numerical models will be validated with the ground-based pressurization and pressure 
control test results. Detailed microgravity simulation case studies for various mission applications can be performed 
with the validated numerical models. The simulation results can then be used to extrapolate and optimize the 
ground-tested spray-bar TVS design for successful microgravity operations. 

To this end, in this paper, we will describe the development of a comprehensive two -phase CFD model for 
simulating the pressurization and pressure control of the Multipurpose Hydrogen Test Bed (MHTB) at NASA 
Marshall Space Flight Center (MSFC), as well as its validation by MHTB ground-based testing results. The 18-m 3 
(693 fi 3 ) storage tank employed in MHTB is representative in both size and shape of a fully integrated space 
transportation vehicle liquid hydrogen (LH2) propellant tank. The tank was designed to accommodate various 
components associated with different CFM concepts including a spray-bar TVS systenr' for active pressure control 
and a multilayer Insulation (MLI) system for passive thermal isolation 2 . 

Numerous numerical models with different degrees of sophistication have been developed to study storage 
tank operations in lg and/or microgravity 4 . However, many of these models have only studied tank self- 
pressurization " 7 . Several have investigated tank pressure control, but mainly for applications where an axial forced 
jet is used to perform thermal destratification through liquid mixing 8 " 11 . In MHTB, storage tank pressure control is 
accomplished through the action of the longitudinal spray-bar TVS that provides simultaneous mixing and cooling 
in the liquid and droplet spraying in the ullage region. Comprehensive models of gas-droplet interaction have been 
developed for thermal management of electronic equipment, evaporative mist flow heat exchangers 12 , and 
combustion engines 13 " 16 . However, applications of such models to cryogenic storage tanks have been quite scarce 
and are limited to lumped models developed by NASA Marshall’ and Rockwell Aerospace 17 and a Flow3D CFD 
model developed through a NASA-Boeing collaboration 18 . The latter predicted the experimental pressurization and 
pressure reduction trends well, but with an average 26% under-prediction of experimental pressure rise and 50% 
under-prediction of experimental pressure drop rates. 

In the two-phase CFD model of the MHTB test tank presented in this paper, the Volume -of-Fluid (VOF) 
technique is used to capture the shape and evolution of the interface. A kinetic-based mass transfer submodel is 
employed to evaluate evaporative -condensing interfacial mass transfer. An Euler-Lagrange approach is utilized to 
track the spray droplets and capture the interaction between droplets and the continuous phase (vapor). By coupling 
the droplet and VOF models, it is possible to perform particle tracking in the vapor, removing particles from the 
vapor domain when they reach the interface and then adding their contributions to the liquid. To account for changes 
in tank pressure during the spraying cycle, an in-house droplet-ullage heat and mass transfer model is also 
developed. The kinetic-based Schrage equation is used to calculate the evaporative/condensing mass transfer from/to 
the droplet. Droplet temperature and size changes are calculated and tracked during the simulations and appropriate 
source terms are included in the governing equations to account for droplet-ullage mass, momentum, and energy 
transfer due to phase change. 

In what follows, the important aspects of the mathematical model will be described first, and then the flow, 
temperature, and interfacial mass flux predicted by the model will be presented. Case studies are performed for the 
50% and 90% fill levels and the predicted ullage pressure and temperatures are compared against their experimental 
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MHTB counterpart for model validation. Simulation results for droplet/ullage heat transfer only and droplet/ullage 
with both heat and mass transfer are compared to underscore the important contribution of droplet 
vaporization/condensation to the evolution of tank pressure during the pressure control interval. 

II. Mathematical Model 


A. Governing Equations 

The geometry and computational grid for the MHTB cryogenic storage tank partially filled with liquid 
hydrogen are shown in Fig. 2. This is a 2D axisymmetric grid that was used in the tank self-pressurization study. 
The grid is shown for a 50% tank fill ratio. A 3D 90° sector grid, created from the 2D grid by revolving it along the 
tank centerline, was used to simulate spray cooling of the MHTB tank. Fluid flow and heat transfer in the tank are 
described in terms of the continuity, Navier-Stokes and energy equations for both phases: 

dp 

-^ + V(pv) = 0, (1) 


dt 


— (pv) + V(pvv) = -Vp + V[p e// (Vv + Vv r )] + pg + F, 


vol> 


- ( pE ) + V(v(p£ + p)) = V(k e// vr) + s h . 


(2) 

(3) 


In the present study, the liquid phase is treated as incompressible with variable temperature -dependent 
properties, except for the density. The liquid density is allowed to vary linearly with temperature in the body force 
term of the momentum equation according to the Boussinesq approximation. The vapor is modeled as a 
compressible ideal gas. All the thermophysical and thermodynamic properties of the fluids are taken from the NIST 
Chemistry WebBook 20 at saturation conditions. 

In this study, the movement of the interface is captured diffusely using the Volume of Fluid (VOF) method, 
as promogulated by Hirt and Nichols 21 . Interfacial energy, momentum and mass balances are applied using source 
terms in the diffuse interfacial region. 


B. VOF Model 

In the VOF method, a volume fraction is defined in each cell such that the volume fractions of all the 
phases sum to unity. In the cell, the change in the interface can be tracked by solving a continuity equation for the 
volume fraction of the q h phase: 


Pq 


— (a q p q ) + V-(a q p q v q ) = S aq 


where the volume fraction for the primary phase is determined from: 


I 

9 = 1 


a q = 1. 


(4) 


(5) 


In the VOF method, the field variables and properties are defined in terms of the volume fraction, which for 
a general system can be written as: 


= 1' 
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In this fashion, the continuity, momentum and energy equations, as described by Eq. (1) - (3), can be 
solved throughout the domain for the temperatures and velocities in the two phases. In the VOF model, energy (E) 
and temperature (7) are treated as mass-averaged variables: 

_ Zg = l a qPqEq , , 

E?=l «,P, ' 

where E q is based on the specific heat of the q ,h phase and the shared temperature. 

Evaporation and condensation at the interface are modeled as a source term in the continuity equation for 
the volume fraction (Eq. 4), i.e.: 

S aq =m r A i , (8) 

where A. is an interfacial area density vector, and m ( is a mass flux vector, which for near equilibrium conditions 
can be determined based on the Schrage 22 equation: 
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Here a is the evaporation efficiency (a value of 0.001 was used in this study); M is the molar mass of the fluid; R is 
the universal gas constant; P, and l\ are, respectively, the interfacial and vapor pressures (it was assumed that P, = 
P sat ); Tj and T v are, respectively, the interfacial and vapor temperatures (it was assumed that T t = T v = T sat at the 
interface). Finally, A, is defined as: 


Aj = |Va|, (10) 

where a is the volume fraction of the primary phase. 

In the present implementation, the surface tension forces at the interface are modeled via the Continuum 
Surface Force (CSF) model of Brackbill et al. 23 . In this model, the surface tension forces at the interface are 
transformed into a volume force ( F vo/ ), which is added as a source to the momentum equation: 


v - ' aiPitiiVaj 

= L aii T7~ 


cCiPihjVoCj + ajpjhjVa i 


pairs ij,i<j 2 (Pi + Pj) 

where h t is the surface curvature calculated from the local gradients in the surface normal at the interface: 

hi = V ■ n. 


( 11 ) 


(12) 


C. Turbulence modeling 

In this study, the tank turbulence was modeled by utilizing the Shear Stress Transport k-w model of 
Menter 19 (kw-SST). This model is similar to the standard k-co model of Wilcox 24 , but has the ability to account for 
the transport of the principal shear stress in adverse pressure gradient boundary layers. The model is based on the 
assumption of Bradshaw et al. 25 that the principal shear stress is proportional to the turbulent kinetic energy, which 
is introduced into the definition of the eddy-viscosity. These features make the kw-SST model more accurate and 
reliable for a wider class of flows than the standard k-co model 26 . 

In the VOF model, continuity of the turbulent quantities is inherently assumed since one set of equations 
for the turbulent kinetic energy and dissipation rate is solved for both phases throughout the domain, with properties 
varying according to the local volume fraction value. 


D. Lagrangian Spray model 

A customized Lagrangian Spray model of the ANSYS Fluent CFD code was utilized for simulating cooling 
of the MHTB tank by spraying cold liquid in a vapor region. This model uses the Euler-Lagrange approach, where 
the fluid phase (ullage) is treated as a continuum by solving the Navier-Stokes equations. The dispersed phase is 
solved by tracking a large number of particles (spray droplets) through the calculated flow field. The droplets 
exchange mass, momentum and energy with the fluid phase. In the Lagrangian spray model, the droplet trajectory is 
calculated by integrating the force balance on the droplet. This force balance equates the droplet inertia with the 
forces acting on the droplet as: 


du.p 

dt 


F D {1 


■ u v ) + 


g(Pp - p ) 
Pp 


+ F, 


(13) 


where U is the fluid phase velocity, Up is the particle velocity, p is the fluid density, p ^ is the density of the 

droplet, F d is the drag force per unit droplet mass, and F is an additional acceleration. Integration of time in 
equation 13 yields the velocity of the droplet at each point along the droplet trajectory. 

The spherical drag law proposed by Morsi and Alexander 27 is applied for the droplet drag force calculation. 
The droplet energy equation is solved in the user subroutine as: 

dT v . . 

mpCv P~dT = h ' A p(. t °° ~ T p) ~ Lrh P’ ( 14 ) 


where m p is the droplet mass, c p is the heat capacity of the droplet, A p is the surface area of the droplet, T m is the 

local temperature of the continuous phase, h is the convective heat transfer coefficient, L is the latent heat, and m p 
is the rate of evaporation/condensation on the droplet, obtained from the Schrage relation (equation 9) in kg/s. 

The temperature of the droplet is obtained from equation 14 using correlation of Ranz and Marshall 2 29 for 
evaluating the heat transfer coefficient: 
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( 15 ) 


Nu = — = 2.0 + 0.6 Re^Pr 1 ' 3 , 

^oo 

where d p is the droplet diameter, k m is the thermal conductivity of the continuous phase, Re d is the Reynolds 
number based on the particle diameter and the relative velocity, and Pr is the Prandtl number of the continuous 
phase. 

The relative Reynolds number is calculated as: 


Re d 


pd p | Up XL 


(16) 


where /. i is the molecular viscosity of the fluid. 

After the droplet energy equation is solved in the user subroutine, the mass, diameter and temperature of 
the droplets are updated. Mass, momentum and energy transfer between the droplets and continuous phase is 
modeled via source terms added to the equations 1-3. Source terms added to the momentum and continuity 
equations for the continuous phase are discussed in AN SYS Fluent theory guide 26 . Source term added to the energy 
equation of the continuous phase is discussed in the results section of this paper. 

The Lagrangian Spray model is coupled with the VOF model via a user subroutine which performs particle 
tracking in the ullage, removes particles from the ullage when they reach the interface, and then adds their 
contributions to the liquid mass, momentum and energy equations through source terms (for those liquid cells near 
the location where the spray drops crossed the interface from vapor to liquid). Because the VOF method produces a 
diffuse interface, the criteria for determining when spray drops crossed into the bulk liquid was defined as the liquid 
drop entering a fluid cell with a liquid volume fraction >0.1. 


III. Numerical Implementation 

The interfacial mass transfer formulation used with the VOF scheme, droplet/ullage heat and mass transfer, 
and the coupling between the VOF and the Lagrangian Spray models have been developed, coded, and implemented 
into a customized in-house version of the ANSYS Fluent version 15.0 CFD code. In the VOF model, the interfacial 
mass transfer is calculated in the user subroutine, and sources are applied to the momentum and relevant scalar 
equations based on the assumption that mass “created” or “destroyed” will have the same momentum and energy of 
the phase from which it was “created” or “destroyed” 26 . All of the two-phase calculations are performed in the user 
subroutine, which is called once per outer iteration before any of the field equations are solved. 

The computational domain for the 50% fill ratio case was discretized using an unstructured mesh of 9,246 
cells, as depicted in Fig. 2. The computational grid for the 90% fill ratio case was created in a similar manner. In 
both the ullage and liquid regions, the conservation equations are evolved in time using a bounded second order time 
stepping routine with a time step size on the order of 5xl0" 2 seconds for self-pressurization cases and lxlO" 4 seconds 
for spray cases. The second order monotone upwinding scheme is used to discretize the convective fluxes in the 
momentum, energy and turbulence equations. The PISO method is used for the pressure-velocity coupling. 

For the volume fraction equation, the Compressive discretization scheme is used. The Compressive scheme 
is a second order reconstruction scheme for the VOF equation based on the slope limiter 26 . In this scheme: 

0/ = 0d+/?Va d , (17) 

where 0/ is the face value of volume fraction, (p d is the donor cell value of volume fraction, /? is the slope limiter, 
and Va d is the donor cell volume fraction gradient. 

In the Lagrangian spray cases, cold liquid enters the liquid phase as a “tiny jet” which is modeled as a 
source term defined in a user subroutine and added to the mass, energy and momentum equations in Fluent. 
Pressure control using spray bar TVS is simulated only for the 90% fill ratio case. In the 90% fill ratio case there are 
38 liquid jets in the liquid region and 4 spray injections in the vapor region. Both liquid jets and spray injections 
have the same uniform temperature of 20.85 K and a liquid flow rate of 6.047 kg/s. The spray injection type is a 
plain orifice atomizer with four particle streams per injection. Two-way coupling between the droplets and 
continuous phase is enabled. In order to model heat and mass transfer between the droplets and ullage, first a value 
for the mass transfer rate is calculated in the subroutine. Then the droplet energy equation is solved, and the droplet 
temperature is updated. The mass and diameter of the droplet are updated to account for the amount of evaporated or 
condensed mass. Then, corresponding mass and energy sources are added to the ullage. 

Convergence criteria are set to lxlO" 5 * for all equations except the energy equation, for which it is set to 

lxlO" 8 . 
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IV. Results and Discussion 

The case studies presented in this section are based on the MHTB tank self-pressurization and spray bar 
mixing experiments, which were conducted in normal gravity 2 . The MHTB tank consists of a cylindrical mid- 
section with 3.05 m diameter and 3.05 m height and two 2:1 elliptical end caps. Two different tank liquid fill ratios 
(50% and 90%) are considered in the MHTB tank self-pressurization study. In the spray bar mixing study, only the 
90% tank fill ratio is considered. 

The gravitational vector is aligned with the tank’s central axis in the negative axial direction. In the 90% fill 
ratio case, a uniform heat flux of 2.05 W/m 2 was applied at the tank wall to match the average wall heat flux of the 
experiment. In the 50% fill ratio case, a uniform heat flux of 0.89873 W/m 2 was applied at the wall in the vapor 
region, and 2.0841 W/m" in the liquid region. These values are based on previous thermal analysis and multinode 
simulations. Conduction through the tank wall was not considered. The spray bar assembly, which is approximated 
as lying along the tank centerline, is treated as an adiabatic surface. A zero velocity field is used as an initial 
condition for the self-pressurization simulations for both fill ratios. The initial turbulence conditions in the tank are 
set to uniform values of k = lxlO" 6 , m 2 /s 2 ; and to = 100, 1/s. The initial tank temperature is set to a uniform 
saturation value corresponding to the tank pressure in the 90% fill ratio case. In the 50% fill ratio case, an 
experimental tank temperature profile is applied as an initial condition for self-pressurization. The details of the self- 
pressurization cases are presented in Table 1. The temperature and velocity fields at the end of the self- 
pressurization case are transferred from 2D into 3D geometry, and used as initial conditions for the spray bar mixing 
case. The details of the spray bar mixing case are presented in Table 2. The initial tank pressure matched the 
experimental one for each case as given in Tables 1 and 2. 


Table 1: Experimental conditions for tank self-pressurization tests 


Case # 

Initial liquid fill, % 

Initial pressure, kPa 

Final Pressure, kPa 

Test Duration, hours 

1 

50 

111.6 

137.9 

13.9 

2 

90 

113.6 

138.4 

6.9 

Table 2: Experimental conditions for spray-bar mixing test 

Case # 

Initial liquid fill, % 

Initial pressure, kPa 

Final Pressure, kPa 

Test Duration, s 

1 

90 

138.4 

131.4 

70 


A grid comparison study was performed for the 50% fill ratio self-pressurization case. The results for the 
tank pressure evolution from the three tested grids are presented and compared with experimental data in Fig. 1. The 
results for the interfacial mass transfer rate from three different grids are compared in Fig. 1 as well, but no 
experimental data on the interfacial mass transfer was available for comparison. All three grids have grid refinement 
near the wall and in the interfacial area. The coarse grid consists of 6,359 cells. The medium grid has 9,246 cells. It 
is similar to the coarse grid, but has more refinement in the interfacial area. It is shown in Fig. 2. The fine grid, 
which consists of 13,482 cells, is further refined at the interface as compared to the medium grid. All three grids 
predict very similar tank pressures and interfacial mass transfer rates. The fine grid simulation required 96 hours 
using 8 AMD Opteron™ 6100 Series (“Magny-Cours”) processors to complete 5000 seconds of self-pressurization. 
The medium grid simulation had approximately ten times shorter running time, compared with the fine grid. The 
coarse grid simulation was approximately twice as fast as the medium grid one. Simulations with the medium and 
coarse grids were run on the same number of processors as the fine grid simulation. Since the results of the self- 
pressurization study will be applied as initial conditions in the spray bar mixing study, where more refinement at the 
interface is beneficial, the medium grid, rather than the coarse one, was selected for this study. A grid with a similar 
level of refinement at the tank wall and at the interface was created and used for the 90% fill ratio self-pressurization 
cases. 

Laminar and turbulent approaches were compared for modeling self-pressurization of the MHTB tank with 
a 50% tank fill ratio. In the turbulent case, the Shear Stress Transport k-co model of Menter 19 was used. The tank 
pressure and interfacial mass transfer rate evolutions from these two models are compared in Fig. 3. The 
experimental tank pressure is also plotted in Fig. 3, but no experimental data on the interfacial mass transfer was 
available for comparison. The tank pressure predicted by the laminar VOF model is in excellent agreement with the 
experiment. This model predicts some evaporation (positive values of the mass transfer rate) at the interface for the 
first 1000 seconds of the tank self-pressurization period. It predicts condensation at the interface for the rest of the 
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simulation time. The turbulent VOF model significantly under-predicts the experimental pressure, and predicts 
evaporation at the interface for the whole duration of tank self-pressurization. Temperatures along the vertical rake, 
located at 0.74168 m from the tank centerline, are compared at the end of the self-pressurization period from the 
laminar and turbulent VOF models and from the experiment in Fig. 4. The turbulent VOF model predicts much 
lower temperatures in the vapor and interfacial regions at the end of self-pressurization, compared to the laminar 
VOF model and experimental data. The laminar VOF model, although it somewhat over-predicts the experimental 
vapor temperatures, agrees with the experimental slope of the temperature distribution in the vapor and with the 
experimental temperatures near the interface. Temperatures and streamlines in the tank from the laminar and 
turbulent VOF models at the end of self-pressurization are compared in Fig. 5. The laminar VOF model predicts 
higher temperatures in the vapor and less mixing in both vapor and liquid at the end of the self-pressurization period, 
compared to the turbulent VOF model. 

The laminar VOF model predicts higher tank pressure and vapor temperatures than the turbulent VOF 
model, despite the fact that the laminar model predicts vapor condensation at the interface during most of the tank 
self-pressurization period. Turbulent VOF model predicts evaporation throughout the self-pressurization period. 
This can be explained by the fact that the laminar VOF model predicts that less heat is transferred in the vapor 
region to the interface due to reduced mixing in the vapor. Thus, less heat is also transferred through the interface 
into the liquid region. The laminar VOF model predicts that more of the heat imparted by the wall stays in the vapor 
region, resulting in higher tank pressure and temperatures. In the turbulent VOF model, more heat is transferred to 
the interface due to increased turbulent mixing in the vapor. Thus, more heat also goes through the interface into the 
liquid. This leaves less heat in the vapor region and results in a lower tank pressure and vapor temperatures. In other 
words, in the turbulent VOF model, more heat is transferred through the interface into the liquid because the 
effective thermal conductivity (molecular plus turbulent) in this region is several orders of magnitude higher than its 
molecular counterpart in the laminar VOF case. This is shown in Fig. 6, where the effective thermal conductivity 
from the turbulent VOF case and the molecular thermal conductivity from the laminar VOF case are plotted along 
the interface at the end of tank self-pressurization. 

Despite the fact that the natural convection flow inside the tank is in the turbulent regime (liquid Ra = 
6.32xl0 14 ; vapor Ra = 1.83 xlO 13 ), the laminar VOF model seems to be better suited for modeling self-pressurization 
of the cryogenic storage tank in normal gravity. The turbulent VOF model implemented in the ANSYS Fluent CFD 
code under-predicts the MFITB tank pressure and vapor temperatures. This may be due to the fact that the turbulent 
damping effect of the interface is neglected in the turbulent VOF model, where continuity of turbulent quantities is 
assumed. 

The MFITB tank self-pressurization study using the laminar and turbulent VOF models was repeated for the 
90% tank fill ratio case. The tank pressure and interfacial mass transfer rate evolutions from these two models are 
compared in Fig. 7. The experimental tank pressure is also plotted in Fig. 7, but no experimental data on the 
interfacial mass transfer was available for comparison. Similar to the 50% tank fill ratio cases, the turbulent VOF 
model significantly under-predicts the experimental pressure. The turbulent VOF model predicts initial condensation 
(negative values of the mass transfer rate) and evaporation at the interface after 8000 seconds of self-pressurization. 
The results of the laminar VOF model agree well with experimental pressures for the first 7000 seconds of self- 
pressurization. After 7000 seconds, and until the end of self-pressurization, the laminar VOF model somewhat over- 
predicts the experimental pressures. This model, just like in the 50% fill ratio case, predicts a short period of initial 
evaporation followed by condensation at the interface for most of the self-pressurization period. The difference in 
the laminar VOF model performance for the 90% fill ratio, compared to the 50% fill ratio case, can be explained by 
the fact that in the 90% fill ratio case, turbulence in the bulk liquid (which is not accounted for in the laminar VOF 
model) affects heat transfer in the tank more. In the 50% fill ratio case, heat transfer in the tank is less affected by 
the turbulence in the bulk liquid because of the constraint imposed by the smaller volume of liquid. 

The MFITB tank spray bar mixing simulation was performed for the 90% tank fill ratio case. The Shear 
Stress Transport k-co turbulence model of Menter 19 combined with the VOF model was used in this case. The tank 
pressure evolution from CFD is compared with experimental data in Fig. 8. The interfacial mass transfer rate 
evolution predicted in CFD is presented in Fig. 8, as well; however, no experimental data was available for 
comparison. The CFD spray model presented in Fig. 8 only includes heat transfer between droplets and ullage and 
assumes no droplet-ullage mass transfer. With this limitation, the CFD spray model predicts no pressure reduction in 
the tank for the first 50 seconds of spray cooling. After 50 seconds of spray cooling, this model predicts that the tank 
pressure begins to drop, but at a much slower rate than in the experiment. The CFD model predicts an initial 
increase in the evaporation rate at the interface, and reduced evaporation after the pressure in the tank starts to 
decrease. 
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The CFD spray model without droplet-ullage mass transfer is not capable of predicting the correct rate of 
tank pressure reduction during spray cooling. Therefore it was necessary to develop and implement a model that 
includes both heat and mass transfer between droplets and ullage. The implementation and results of this model are 
discussed below. The Schrage relation (Eq. 9) was utilized to calculate the mass transfer rate between the droplet 
and ullage. Details of this model are presented in the Mathematical model section of this paper. As a first step in 
model implementation, a test study was performed in order to identify the correct mass and energy sources that are 
needed to couple the droplets and continuous phase (ullage) during mass transfer. In the test study, a simplified 
closed volume filled with compressible hydrogen vapor was considered. Two sides of the closed volume are square 
walls with 9 by 9 cells of orthogonal mesh. The extent of the domain is 10 times its height. The initial vapor 
temperature was set to a uniform value of 22 K. In order to test the source terms added to the continuous phase 
during mass transfer, first two liquid hydrogen droplets were released at the center of the square wall and the 
simulation was run for 10 time steps. The initial temperature of the droplets at injection was set to 21 K. The 
resulting vapor and droplet temperatures, vapor mass change and vapor energy sources are presented in Figures 9, 
10 and 11, respectively. As shown in Fig. 9, the droplets only travelled through 2 computational cells during the 
simulation. After 10 time steps, the vapor temperature decreased in those cells that the droplets travelled through. 
This is due to the heat transfer between the hotter vapor and cooler droplets. The temperature of the droplets 
increased, as expected. Figure 10 shows that the vapor mass (monitored in Fluent) decreased during the simulation, 
which indicates vapor condensation on the droplets. The vapor mass source due to droplet/vapor mass transfer is 
implemented correctly if the vapor mass change matches the total amount of mass transfer from the droplets 
calculated in the user subroutine. This was tested, and the results are shown in Fig. 10. Here, the vapor mass change 
monitored in Fluent is compared with the total droplet mass change (with a negative sign) calculated in the user 
subroutine. The fact that the two results match indicates that the vapor mass source due to droplet/vapor mass 
transfer is implemented correctly. 

According to the ANSYS Fluent Documentation 26 , the heat transfer from the continuous phase (vapor) to 
the discrete phase (droplet) is computed by examining the change in thermal energy of a droplet as it passes through 
each control volume, as: 
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where m p 0 is the initial mass flow rate of the particle injection (kg/s); m p 0 is the initial mass of the particle (kg); 
m p . is the mass of the particle on cell entry (kg); m p ( is the mass of the particle on cell exit (kg); c p is the heat 

capacity of the particle (J/(kg-K)); T p .^ is the temperature of the particle on cell entry (K); T p is the temperature 
of the particle on cell exit (K); T re f is the reference temperature for the enthalpy equation (K); L re f is the latent heat 
at reference conditions (J/kg). 

The first part on the right hand side of Eq. 18 (in square brackets) is the correction for the particle mass 
source being added at a temperature other than the Fluent reference temperature. For simplicity, we will call this 
source “correction for the reference temperature”. The second part on the right hand side of equation 18 is the 
energy source representing the latent heat addition. This source has to be added at the reference temperature. It is not 
clearly stated in ANSYS Fluent documentation what energy sources are added to the continuous phase internally if a 
mass source is added in a user subroutine. In order to investigate this issue, the sources presented in Eq. 18 are 
plotted in Fig. 11 for the duration of the spray test case. It is clearly seen that when only the latent heat addition 
source is added in the subroutine, the total source added in Fluent also includes the correction for the reference 
temperature part; therefore, there is no need to add this part in the subroutine. 

To make sure that this model is also applicable for multiple droplets injected every time step and that it can 
be run for a longer time, the spray test case was rerun with two droplets injected every time step for 100 time steps. 
Vapor and droplet temperatures predicted in this case are presented in Fig. 12. Similar to the results of the previous 
test case shown in Fig. 9, here vapor is being cooled by the droplets and the droplets are being heated by the hotter 
vapor, as expected. Based on the analysis of the spray test cases presented above, the implementation of the 
droplet/ullage heat and mass transfer model is verified and can be applied to simulate the MHTB spray bar mixing 
experiment. 
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The results of the MHTB spray bar mixing simulation with the droplet/ullage heat and mass transfer model 
are presented in Fig. 13. Here, the tank pressure evolutions are compared between the two models, with and without 
droplet/ullage mass transfer, and the model results are also compared with experimental data. The interfacial mass 
transfer rate evolutions predicted by the two models are compared in Fig. 13, as well; however, no experimental data 
was available for comparison. The model with droplet/ullage heat and mass transfer predicts a tank pressure drop 
during spray cooling that compares well with experimental data for the first 25 seconds of the simulation. After 25 
seconds of spray, this model somewhat over-predicts the rate of the pressure drop in the tank. The model with only 
heat transfer between the droplets and ullage predicts no pressure reduction in the tank for the first 50 seconds of 
spray, which disagrees with the experiment and the model with droplet/ullage heat and mass transfer. The model 
with heat and mass transfer between droplets and ullage predicts less evaporation at the interface compared with the 
model that only includes heat transfer. This is because, in the model with only heat transfer modeled between the 
droplets and ullage, the vapor is not cooled by evaporating droplets, and hotter vapor at the interface results in more 
evaporation than in the case with droplet/ullage heat and mass transfer. This is illustrated in Fig. 14, where 
temperature and velocity fields at the mid-plane section of the tank predicted by the two models are shown after 50 
seconds of spray. In the experiment the spray is on for 70 seconds, but significant tank pressure reduction is already 
achieved at 50 seconds of spray. Droplet distributions from the two models colored by droplet temperatures are 
shown in Fig. 14, as well. Higher droplet temperatures and higher temperatures in the ullage are observed in the case 
with no droplet/ullage mass transfer. In the case with droplet/ullage mass transfer, evaporating droplets cool down 
the ullage. This case results in lower droplet and ullage temperatures than the case without droplet/ullage mass 
transfer. Evaporating droplets lose their mass and momentum, and do not travel through the vapor as far as the 
droplets that do not evaporate, as can be seen from the velocity vectors shown in Fig. 14. 

A droplet/ullage heat and mass transfer model was developed and implemented with Lagrangian spray and 
the VOF models of ANSYS Fluent commercial CFD code. This model was utilized for simulating the MHTB spray 
bar mixing experiment. The model with droplet/ullage mass transfer predicted a pressure drop in the tank that 
matched the experimental one. The CFD spray model without droplet/ullage mass transfer predicted a delayed 
pressure reduction in the tank and a slower pressure reduction afterward, compared to the experimental data. This is 
due to a larger amount of heat transferred by the droplets to the interface, which resulted in a high evaporation rate . 
It is necessary to include droplet-ullage mass transfer into a CFD model in order to accurately predict the spray 
cooling of cryogenic storage tanks. 


V. Conclusion 

A two-phase CFD compressible VOF model for the self-pressurization and pressure control inside the 
cryogenic storage tanks was presented. Two different tank fill ratios, 50% and 90%, were considered. The laminar 
and turbulent approaches to modeling the liquid, ullage, and interfacial regions were compared. The model was 
validated against the experimental data provided from the MHTB self-pressurization and spray bar mixing 
experiments. Simulation results for droplet/ullage heat transfer only and droplet/ullage with both heat and mass 
transfer were presented and compared to underscore the important contribution of droplet vaporization/condensation 
to the evolution of tank pressure during the pressure control interval. 

In the case of the tank self-pressurization the laminar VOF model resulted in tank pressures and 
temperatures that agree well with the experimental ones for the whole duration of simulation in the 50% tank fill 
ratio case. In the 90% tank fill ratio case, the tank pressure predicted by the laminar VOF model agreed well with the 
experiment for the first 7000 seconds of self-pressurization. This model predicted higher than experimental tank 
pressure for the rest of the self-pressurization period. The turbulent VOF model significantly under -predicted the 
tank pressure and temperatures in the ullage compared to the experimental data. This is due to excessive mixing in 
the ullage and over-prediction of the amount of heat transfer through the interface produced by this model. 

A CFD model for cooling of cryogenic storage tanks by spraying cold liquid in the vapor region was 
developed, and its results were presented. The Euler-Lagrange approach implemented in the spray model of the 
ANSYS Fluent CFD code was utilized for tracking the spray droplets, and for modeling the interaction between the 
droplets and the continuous phase (ullage). The spray model was coupled with the VOF model via a user subroutine 
which performs particle tracking in the vapor, removes particles from the vapor domain when they reach the 
interface, and then adds their contributions to the liquid through source terms. Two different versions of this model 
were compared. The first version included only heat transfer between the droplets and ullage. This version is part of 
the ANSYS Fluent code. In this version, a delay was predicted in the tank pressure reduction during spray cooling as 
compared to the experiment for the first 50 seconds of spray, and the tank pressure decreased slower than in 
experiment for the rest of the simulation. This is due to a large amount of heat transferred by the droplets to the 
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interface, which resulted in a higher evaporation rate. In the second version, a model for mass transfer interaction 
between droplet and ullage due to droplet phase change was developed and incorporated into Fluent via a user 
subroutine. The model with droplet/ullage heat and mass transfer resulted in a more accurate prediction of the tank 
pressure decrease during spray cooling of the cryogenic storage tank. This model also predicted a lower amount of 
heat transferred by the droplets to the interface, which resulted in a lower evaporation rate than predicted by the 
model with only droplet/ullage heat transfer. 

The results of this study show that it is important to include the droplet/ullage heat and mass transfer when 
simulating active cooling of cryogenic storage tanks using spray. 
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Figure 1: Pressure and interfacial mass transfer rate (positive for evaporation) evolutions during self-pressurization 

of MHTB tank (50% fill ratio): Grid Comparison 
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Figure 3: Pressure and interfacial mass transfer rate (positive values represent evaporation) evolutions in the tank 
during self-pressurization of MHTB tank (50% fill ratio): Effect of turbulence modeling with the VOF model 
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Figure 4: Comparison between laminar and turbulent VOF cases (50% fill ratio): Temperatures along the vertical 
rake (0.74168 m off the tank centerline) at the end of tank self-pressurization 
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Figure 5: Comparison between laminar and turbulent VOF cases (50% fill ratio): Tank temperature field and 

streamlines at the end of tank self-pressurization 
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plotted along the interface at the end of tank self-pressurization 
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Figure 7: Pressure and interfacial mass transfer rate (positive in case of evaporation) evolutions in the tank during 
self-pressurization of MHTB tank (90% fill ratio): Effect of turbulence modeling with the VOF model 
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Figure 8: Pressure and interfacial mass transfer rate (positive in case of evaporation) evolutions in the tank during 

spray cooling of MHTB tank (90% fill ratio) 
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Figure 9: Cell by cell analysis of the droplet/ullage condensation test case ran for 10 time steps with only two 

droplets injected at the first time step (At = 0.0001s) 



Figure 10: Comparison between vapor mass change in Fluent and total droplet mass change in the subroutine (sign 

changed) during 10 time steps 
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Figure 11: Analysis of the enthalpy source added to the vapor due to condensation on the droplets during 10 time 

steps 
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Figure 12: Cell by cell analysis of the droplet/ullage condensation test case ran for 100 time steps with two droplets 

injected every time step (At = 0.0001s) 
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Figure 13: Pressure and interfacial mass transfer rate (positive in case of evaporation) evolutions in the tank during 
spray cooling of MHTB tank (90% fill ratio) - results of the model with droplet/ullage mass transfer added 


18 

American Institute of Aeronautics and Astronautics 


mass transfer rate at the interface, kg/s 


T_droplet, K 





V, m/s 

1.00 

0.90 

0.80 

0.70 

0.60 

0.50 

0.40 

0.30 


0.20 


0.00 


Figure 14: Vapor and droplet temperatures and flow field in the MHTB tank after 50 seconds of spray cooling (90% 
fill ratio) - comparison between the models with and without droplet/ullage mass transfer 
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